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Abstract 

This paper presents an analysis procedure for experimental data using theoretical 
functions generated by Monte Carlo. Applying the classical chi-square fitting pro- 
cedure for some multiparameter systems is extremely difficult due to a lack of an 
analytical expression for the theoretical functions describing the system. The pro- 
posed algorithm is based on the least square method using a grid of Monte Carlo 
generated functions each corresponding to definite values of the minimization pa- 
rameters. It is used for the E742 experiment (TRIUMF, Vancouver, Canada) data 
analysis with the aim to extract muonic atom scattering parameters on solid hydro- 
gen. 
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1 Introduction 



For a wide range of physical problems the only applicable way to compare 
experiment with theory is via the Monte-Carlo (MC) method. Problems of that 
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type are often multiparameter, with nontrivial interdependcncies between the 
parameters, such as averaging arising from spatially discrete effects. Thus, only 
an exact simulation of the experimental system allows us a possible analysis, 
and thus we rely on the MC method. 

However, the MC method has several limitations, mostly related to calculation 
time and the nature of random sampling. Provided the simulation has been 
correctly established to analyze all competitive processes, the modeling of 
events which are fairly rare and hidden inside many other processes often 
requires long calculation time to establish sufficient statistics for comparison. 
The intrinsic nature of random numbers, and the generators presently in use, 
give that any two simulations of the same system can give different results. 
Normally, for large numbers of generated statistics those differences are small 
and completely insignificant. Such limitations are not very important when we 
make a single simulation for one experiment, i.e., when there are no variable 
parameters and one output result is sufficient. However, the limitations are 
amplified when we apply the MC method to a fitting procedures. 

As is fairly well known, finding the best fit parameters describing a multitude 
of data requires repeated calculation of some statistical estimator. Most often, 
the estimator is which is defined as the difference between the theoretical 
description and the experimental data, see Eq. (3). The theoretical description 
of the data depends on several parameters and thus is calculated as a 
function of those parameters. Finally, the result is the set of parameters for 
which the is minimal. 

Classical fitting algorithms, e.g., the minimization package MiNUlT [?], cal- 
culate from the model parameters (see Fig. 1). When the minimum 
depends on two or more variables, the error determination on the parameters 
as well as the study of the possible occurrence from several minima require 
calculations of several thousands theoretical functions, and hence, the cal- 
culations becomes extremely time-consuming. However, the MC evaluation 
of the theoretical function, just for one set of parameters, is very time ex- 
haustive (measured in hours or even days): thousands of iterations are not 
possible. Even if the calculation time were acceptable, the intrinsic nature 
of MC simulations makes such an approach impossible since instabilities will 
arise resulting from the statistical nature of the results. ^ . 

When fitting, the minimization procedure examines the behaviour of differ- 
ences in for differing values of the parameter set. The minimization pro- 
cedure then calculates the internal gradient of x^ and uses it to control the 

^ For example, MiNUlT examines whether the theoretical function is time- 
independent. The theoretical function for a given parameter set is evaluated twice. 
If the resulting values are different, the theoretical function is qualified as a time- 
dependent and the minimization procedure is suspended. 
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minimum searching procedure. The gradient is obtained from a set of partial 
derivatives for each variable parameter where the derivatives are calculated nu- 
merically from difference quotients. Normally, the minimum should be reached 
when all the gradients converge to zero. Clearly, the statistical fluctuations of 
the MC method can cause entirely false gradients, and thus such a minimiza- 
tion procedure is not suitable for our problem. 

The work of Zech [?] presents methods for comparing MC generated his- 
tograms to experimental data when the analytic distribution is known. In 
our case we compare the experimental data with MC simulations [?] which 
including all parameters of the apparatus, such as the spatial separation of 
processes, detector resolution, dead time, etc, and therefore, we can directly 
compare experimental and Monte Carlo spectra. 

If wc ask "given a data distribution, and a set of MC distributions, what is 
the best estimate of the fraction of each MC distribution present in the data 
distribution?" a standard set of subroutines [?], are available to solve that 
problem. However, the experimental histograms in our case are not equivalent 
to a summing of discrete MC spectra, and the method above cannot be ap- 
plied. An approach [?] similar to ours was used to determine the muon energy 
distributions following muonic capture and atomic cascade using time of flight 
methods, although no detailed description of the method is available. The aim 
our paper is to describe such algorithms and to show via example that it is 
fully applicable. As an example the scattering of muonic atoms on a structure 
of crystalline hydrogen is presented. 



2 Description of the Method 

2. 1 Modified fitting procedure 

We proposed a modiflcation in the calculation of the theoretical functions 
M{'p) which describe the data for a given parameter vector p — {pi,p2, ■ ■ ■ ,Pn)- 
Before fltting, one generates a set of theoretical functions {M} = {M^^s,...,(j), 
My^gr (^/, . . .} for all permutations of a chosen discrete set of parameter val- 
ues (^1,^2! • • • )Pn) where j,5, . . . ,(f) are the function indexes in the set. The 
parameter vector p is allowed to assume only discrete values which gives the 
grid of theoretical functions a size j x S x . . . x (f), and means that the time- 
consuming calculations are only executed for a select and limited parameter 
set {p} = {{Pi,pI, ■ ■ ■ ,pt)i (Pi 1P2 5 • • • )i • • •}• The resulting set of function 
values, {M}, is used to calculate any theoretical function M for any arbitrary 
parameter vector p (provided all Pi values in p are between some calculated 
values of p] and p] contained in the grid) using an interpolation procedure 
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Fig. 1. Scheme of the fitting procedure using grids. The modifications are marked 
with a thick line. The dashed line denotes the procedure defined by the users. 

described in Sec. 2.3. Once the {M} set is known, the interpolation is rela- 
tively fast, and the results, M{p), can be used to calculate the x^- Note that 
none of the above precludes M from depending on other variables, such as 
time or space, and hence the generated M^^s,...,4> could just as well be written 
M^,s,...,<j>{t, x), so the generated functions may very well, themselves, be multi- 
dimensional. Figure 1 presents schematically the fitting procedure with these 
modifications. 

The number of functions in the set {M} depends on each analysis case and 
should depend on the behaviour of the function M{p) for a given parameter pi. 
One should note that using too few grid points will give only a weak expression 
of the function's behaviour, whereas using too fine a division will, for small 
parameter changes, falsify the gradient calculations due to the statistical MC 
fluctuations. Properly chosen distances between grid points eliminates the 
statistical fluctuations of the theoretical function because the values between 
grid points are interpolated. 



2.2 Description of the calculation 



We choose MC statistics on average about ten times greater than the statisti- 
cal uncertainty in experimental data (less than that and we are insensitive to 
our parameters while fitting; more and we use more MC time for essentially 
no gain in sensitivity). Therefore, we neglect the statistical errors connected 
with the MC and use the classical definition where the fits of the analytical 
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ftmctions arc applied. Very sophisticated definitions of x^, including MC sta- 
tistical fluctuations, are presented in Ref. [?], however, they are most useful 
in the case where experiment and simulation have similar statistics. 

It is possible to use many sets of data from different experimental conditions, 
provided they can all be modeled by the same p parameters, and define the 
total as the sum of the individual xi calculated separately for a single data 
set k. Thus, the total x^, calculated when we perform simultaneously fits to 
m sets of data, is: 

X^^m^^Wfe-Xfe (1) 

k 

where m is the number of fitted histograms, k the histogram index running 
over a single set of data, and Wk the corresponding weight. The factor m is 
used to compare the number of degree of freedom since the chi-squares are 
non- normalized. 

The weights, w^, are calculated as a count ratio in each histogram relative to 
the total counts in all histograms, such that histograms with more counts give 
greater share in the total x^: 

= ^ (2) 

where A*"* is the number of events in channel i of the experimental spectrum. 
The partial x^ is calculated as: 

, [c, ■ M\k) - N\k)f 

where Ck is a factor matching the k^^ experimental TV* with its corresponding 
MC histograms M* and is given by: 

Y.N\k) 



2.3 The interpolation method 



The grids are generated only for a finite and discretized set of parameters for 
all permutations of the parameters. However, as follows from the minimiza- 
tion procedure, theoretical functions are necessary from a continuous param- 
eter space (pi,P2, • • • ,Pn) and the interpolation procedure using the grids is 
applied to generate such functions. A visual scheme of the two-dimensional 
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interpolation procedure is presented in Fig. 2, an example taken from Ref. [?]. 
In general, interpolation is only well defined for scalar values. An interpola- 
tion procedure on function can only occurred if we treat it as a set of scalars. 
Therefore the function is given as a table of scalars. For each table value the 
interpolation is performed separately. Then the set of interpolated values give 
the final interpolated function. 

Figure 2 presents the two-dimensional plane for the parameters pi and p2, 
and represents a part of a full grid. The point at (^1,^2) which we wish 
to find the function is marked by the sign 0. The upper index a means that 
the variable may takes values from the continuous spectrum rather than only 
grid point values. The four points with index a, shown with the open circle 
symbol O) ^-^e intermediate points required by the calculation. Around them 
are the four grid points, (pf,P2), (p^^Spf), (pf,P2^^), and (pi^^,P2^^), 
denoted with the filled circle • symbol. They correspond to the grid functions 
Ma,i3, Ma+1,/3, Ma,/3-\-i, and Ma+i,/3-\-i, respectively. The functions M are given 
in value-channel numerical form, i.e., a number of counts for each channel of 
the spectrum. 

The idea of a two-dimensional interpolation relies on first performing a one- 
dimensional interpolations, where a is the grid size in this direction, along the 
directions connecting the grid points to find the function at point (p^,P2)) 

1^*. One dimensional 
interpolations for each p2 



'2 

'c. + l(P§) 



Fig. 2. Outlook view of the two-dimensional interpolations. Symbols are given in 
the text. 
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for example. In this way, values for the intermediate function /ia(P2) de- 
termined for each node a and a + 1. Secondly, a single onc-dimcnsional in- 
terpolation along the horizontal axis (marked with a dashed line) is executed. 
To obtain a temporary interpolating function g{pi), one needs a function M 
in the point (K)P2)- To make the interpolation function M it is necessary to 
repeat the described procedure for each tabulated value of the independent 
variable ^ . The one-dimensional interpolation is described below. 

One uses the following interpolation formula. 



+ C{pl) ■ [K{P2)T + D{p^2) ■ [K+M)T ^ ^ 

and finally: 

(7) 

where A, B, C, and D are the interpolation coefficients. The coefficients A 
and B are defined as 

A{x) = , B{x) " 



Ax Ax 

where x is used as a formal notation for the independent variable, x e {xa, Xa+i) ■ 
In our case x plays the role of the parameters pi or p2. Ax ~ Xa+i — Xa are 
the grid steps. The coefficients C and D are expressed as: 

C{x) = ^(^A^-A)-Ax\ D{x) = ^{B''-B)-Ax^ (9) 



The dependence of the interpolated function or on x in Eqs. (5) or (6) 
is given by a linear dependence on the coefficients A, B, and a cubic depen- 
dence on the coefficients C, D. Thus, this method is called the cubic spline 
interpolation. 

The first step in this method is the calculation and tabulation of the second 
derivative values for all functions {M^'^^} in the grid. If one wants to use Eq. (6) 
the second derivative of h" is required. One obtains the second derivatives by 
solving the following expressions 

M'U, + 4M:,, + M:^,, = 6 • , (10) 



^ The central points and widths of the channels have to be the same for all the grid 
points. 
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This expression is correct only when grid points are spaced equally. Based 
on Eq. (10), a triangular matrix is created and reduced by using a suitable 
numerical algorithm. To solve Eq. (10) one needs boundary conditions. In the 
presented analysis, the second derivatives for the first and last values of M" 
are set to zero, the so-called natural cubic spline interpolation. 

The ability of our interpolation routine to reproduce a function at (a, (3) using 
the four points (ail, /?±1) is shown in Fig. 3 for a case where the two functions 
y(shift, depth) are given. 



7000 1 r 0.2 




2000 -I , , , , , , , , , , , , , , 1 -0,6 

0,5 0,6 0,7 0,8 0,9 1 1,1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 1,9 2 

deep value 



Fig. 3. Graph of two functions, y(shift = —0.1 eV, depth) given by the red 
squares and y (shift = —0.0 eV, depth) with the green squares. The red tri- 
angles show y" (shift = —0.1 eV, depth) , whereas the green ones stand for 
y"(shift = 0.0 eV, depth). The green solid Une is obtained by an interpolation func- 
tion for a shift value of —0.05 eV. 



3 Application of the method 

3. 1 Description of the experiment 

In this section we apply our analysis method to the data obtained in the 
E742 experiment performed at TRIUMF, Vancouver (Canada). The exper- 
iment was dedicated to the study of //-atomic processes occurring in solid 
hydrogen isotopes. It is of particular interest to obtain the characteristics of 
the interacting systems (collision energy of muonic atoms with a crystalline 
structure). We want to reconstruct the energy dependence of the elastic scat- 
tering cross-sections for muonic atoms in the process: dfi + p —>■ dfi + p on 
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1 



crystalline hydrogen at a temperature of 3 K. Theoretical calculations [?] pos- 
tulated that there exists an energy region of abnormally small cross-section 
called the Ramsauer-Townsend (R-T) region. Figure 4 shows this dependence 
and the R-T region is visible. The aim of the measurement was to find exper- 
imentally the R-T region and verify the theoretical prediction. The accuracy 
of the theoretical calculations was low for energies inside the R-T region and 
the R-T minimum was also determined by theoretical calculations of poor 
accuracy. We assumed only that the general shape of the cross-section curve 
was valid. We vary only the depth and the position of the minimum of the 
R-T region. 




2 



0.001 0.01 0.1 1 10 100 

Energy (eV) 

Fig. 4. Elastic scattering cross-sections for djj, + p and R-T region. The essential 
values used in the parametrization are given. 

The most accurate experimental method would be the use of a (selectable) 
monoenergetic beam of muonic d^i atoms, and, by aiming the beam at a thin 
foil of crystalline hydrogen and (like in the Rutherford experiment) detecting 
the intensity and energy of the scattered dji atoms would allow us to determine 
the cross-section as a function of the d/i-atoms energy. The method would 
also give us the scattering angles, and interpreting the data would be easy. 
However, even in this maximally simplified case it is most probable that the 
MC method have to be used. 

Unfortunately, it is impossible to use the direct method because such a source 
of d^jL atoms does not exist and there are no detectors measuring directly the 
energy of muonic d/x atoms. Therefore, in the real experiment we produced 
muonic atoms inside a structure of solid hydrogen, starting a time counter 
on the muon arrival. The energy of the scattered d^ atoms after leaving the 
crystal is measured indirectly via a time of fiight method. The muonic atom 
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flew between two layers placed at some distance between each other in vacuum. 
The first layer is a source of muonic atoms and is treated as an emitter of di- 
atoms. The second layer is covered by neon and is treated as the detector of 
muonic atoms: a d/i atom entering the neon layer transfers the muon to the Ne 
almost immediately yielding an x ray which determined the stop of the time 
counter. Detailed information about the experiment can be found in Ref. [?] 
and references therein. 

To analyze the experiment, we have to use MC simulations because other pro- 
cesses are competing with the scattering, as well as the complications coming 

from the geometry of the experiment. Figure 5 presents the scheme of the 
processes taken into account [?]. The initial time is given when a muon en- 
ters the apparatus. The output of the simulation is a x-ray time spectrum 
(see processes in external layers) which can be directly compared with the 
experimentally measured one. 
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Fig. 5. Scheme of FOW program used to simulate the E742 experiment. Four process 
blocks are shown, namely the processes in external layers, the muon stopping, the 
diffusion, and the muonic processes. 



3.2 The grid construction 



The grid method is illustrated for one chosen set of three experiments with 

different experimental conditions. Since the scattering cross-section does not 
depend on experimental conditions, fitting three different conditions with the 
same parameters should reduce systematics. The theoretical dependence of 
the cross-section from the collision energy was parametrized on two ways. 
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One changes the position of the minimum on the energy axis (varying AE) 
and the depth of R-T minimum (varying s). In our example pi represents the 
shift of R-T energy (AE) as can be seen in Figs. 4 and 6. The parameter p2 
represents a rescahng factor s of the minimal cross-section value (see Figs. 4 
and 8). The function M is the MC time spectrum of muonic atoms reaching 
the neon layer for AE and s. Thus, the grid is defined as {Mab„,s^}- 



3.2.1 Change of the R-T minimum energy 
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Fig. 6. Ai?-parametrization of the R-T position minimum. Curves for some value 
of AE are shown. 

The energy axis is transformed according to E ^ Etr and then a{E) — > a{Etr). 
The energy transformation are done via 



E, 



tr 



E + • AE, 



E- 



for E < E 



R 



■ AE, for E> Er 



Emax—E^ 

Table 1 

Characteristic values used in the cross section parametrization. 
Notation Name Value 



Er Energy of R-T minimum 

Emin Minimum of energy range 
^max Maximum of energy range 



1.7 

0.001 

190.5 



(11) 



Units 
eV 
eV 
eV 



(t{Er) Minimum cross section for the R-T energy 1.13 x 10 cm' 
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where AE is the energy shift, E the unshifted energy, and Eji the original 
vahic of the R-T energy minimum. Limits Emm and E„iax define the range 
where the transformation is apphed. For E = E^ the transformed energy E^^ 
given by Eq. (11) becomes Etr — E + AE, and thus, this parametrization can 
be treated as a shift. Characteristic values for theses variables are given in 
Table 1. 
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Fig. 7. Ma{AE), the time spectra for the shift-parametrization and s 
are shown for the same values as given in Fig. 6. 



1. Curves 



The shift values, AE, were defined as 11 points between —0.5 eV to 0.5 eV, 
in steps of 0.1 eV. These values were chosen as a result of previous tests 
of the experimental data with different shift values. The shape of the cross- 
section curves are presented on Fig. 6. The resulting MC time spectra for such 
parametrized cross-section are presented in Fig. 7. 



3.2.2 Change of the R-T depth minimum 

The curves for the depth parametrization have a rescaled minimum cross- 
section a{Eji) — s X a{Eji). The rescaling took place in the energy range 
0.5 — 6 eV, with 11 values of the rescahng parameter s chosen from 0.5 to 1.1 
by steps of 0.1, with additional values of 1.25, 1.4, 1.7, and 2.0, as can be seen 
in Fig. 8. 

To preserve the smooth form of the cross sections, values within the 0.5- 
6 eV range were not globally scaled by s, but were scaled by a factor which 
ranged from 1 at the borders, to s at the energy of the R-T minimum. The 
s values (except the minimum value) were selected numerically to reproduce 
the characteristic shape of the cross-section function inside the R-T region, 
as shown in Fig. 9. The MC time spectra for the depth parametrized cross- 
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Fig. 8. Depth-parametrization of the R-T minimum depth. Some values of s, which 
concern the minimum cross section a{Eii) are indicated. 



sections are presented in Fig. 9. The functions Ma{^E) and Mq,(s) given 
in Figs. 7 and 9, respectively, were combined and a grid {Mq,^^(A£', s)} was 
created. 




Fig. 9. M^(s), the time spectra for the depth parametrization, without the energy 
shift A£: = 0. 
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4 The results 



An example of the experimental data fits using the MC time spectra (see 
Figs. 7 and 9) is presented in Fig. 10. In this case three sets of experimental 
data (called Expositions 1-3) were fitted. Fits were performed for a number 
of data combinations. 
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Fig. 10. Example of a fit. 

The average values from all possible fit combinations give the final result: 

s = 1.12 ±0.20, A£; = 0.30 ± 0.14 eV. (12) 

The errors are connected with low experimental statistics, some background 
problems, and also with the grid steps (0.10 eV for AE and 0.1 — 0.15 for 
s) and finally with the interpolation procedure. The result means that our 
experimental result confirms the theoretical cross-section energy dependence 
in the R-T region, but indicates that the minimum of the cross-section value 
occurs for energies higher than predicted, at about 2 eV instead of 1.7 eV in 
Fig. 4. The absolute value of the fitted cross-sections agree with the theoretical 
value. 



5 Conclusion 



The method allowed us to perform a correct comparison of experimental data 
with theoretical predictions based on MC calculations. Although the experi- 
mental data was obtained in only a few weeks of muon beam usage, the grid 
construction was a time-consuming step requiring more than six months of 
calculation. The fitting procedure was quick and allowed us to prepare many 
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fits for any combination of the data and perform more complex analysis of the 
data themselves, e.g., x^-contour and error calculations. 

To estabhsh a more precise set of mathematical rules which test the correctness 
of this method, one needs to perform further studies but such was not the aim 
of this work. The grid method, as demonstrated here, is fully acceptable for 
analyzing systems with complex multiparameter dependences. Our procedure 
was internally checked by comparing results of fitting single and summed data. 
The errors of single data fits are bigger but they lie within the range of the 
experimental errors of summed data. 
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